One-way sensitivity analysis - pre-calibration for status quo model

Published

17/04/2023 T15:10

Code
library(here)
library(plotly)
theme_set(theme_minimal())

source(here("02_scripts/01_fun_data.R"))

t_owsa_ini_pop <- read.csv(file = here("01_data/ows_tbl_ini_pop.csv")) %>% select(-X)
ows_tbl_prob <- read.csv(, file = here("01_data/ows_tbl_prob.csv")) %>% select(-X)

1 Initial population

Code
# m_owsa_ini_pop <- matrix(0, 
#                          nrow = nrow(ini_states_tbl_org),
#                          ncol = 8,
#                          dimnames = list(
#                            rows = paste0("ini_pop_",v_state_names),
#                            cols = c("input_low", "input_high",
#                                     "costs_low", "costs_high",
#                                     "deaths_low", "deaths_high",
#                                     "od_deaths_low", "od_deaths_high")
#                          ))
# 
# for (state in v_state_names){
#   
#   t_ini_pop_temp <- ini_states_tbl_org %>% 
#     mutate(basevalue = ifelse(`Variable ` == state, lb, basevalue))
#   
#   mod <- markov_mod(num_cycles = n_cycles,
#            vec_state_names = v_state_names,
#            vec_m_0 = t_ini_pop_temp$basevalue,
#            mat_P = NA,
#            array_P = arr_P,
#            vec_cost_states = v_cost_states,
#            disc_fac = 0.03, t_0_cost = 0, inc = pop_increase_tbl_new$increase)
#   
#   deaths_lb <- mod$m_M[181, "BO_DEATH"]
#   od_deaths_lb <- mod$m_M[181, "BO_OD_DEATH"] + mod$extra_od_deaths
#   costs_lb <- mod$total_net_present_cost
#   
#   t_ini_pop_temp <- ini_states_tbl_org %>% 
#     mutate(basevalue = ifelse(`Variable ` == state, ub, basevalue))
#   
#   mod <- markov_mod(num_cycles = n_cycles,
#            vec_state_names = v_state_names,
#            vec_m_0 = t_ini_pop_temp$basevalue,
#            mat_P = NA,
#            array_P = arr_P,
#            vec_cost_states = v_cost_states,
#            disc_fac = 0.03, t_0_cost = 0, inc = pop_increase_tbl_new$increase)
#   
#   deaths_ub <- mod$m_M[181, "BO_DEATH"]
#   od_deaths_ub <- mod$m_M[181, "BO_OD_DEATH"] + mod$extra_od_deaths
#   costs_ub <- mod$total_net_present_cost
#   
#   m_owsa_ini_pop[paste0("ini_pop_",state), ] <- c(ini_states_tbl_org$lb[ini_states_tbl_org$`Variable ` == state],
#                                                   ini_states_tbl_org$ub[ini_states_tbl_org$`Variable ` == state],
#                                costs_lb, costs_ub, 
#                                deaths_lb, deaths_ub,
#                                od_deaths_lb, od_deaths_ub)
# }
# 
# 
# t_owsa_ini_pop <- as_tibble(m_owsa_ini_pop) %>% 
#   mutate(costs_range = abs(costs_high - costs_low),
#          deaths_range = abs(deaths_high - deaths_low),
#          od_deaths_range = abs(od_deaths_high - od_deaths_low),
#          name = rownames(m_owsa_ini_pop))

1.0.0.1 Costs

Code
t_owsa_ini_pop_costs <- t_owsa_ini_pop %>% 
  select(costs_low, costs_high, name) 
t_owsa_ini_pop_costs$name <- reorder(t_owsa_ini_pop_costs$name, t_owsa_ini_pop$costs_range)

t_owsa_ini_pop_costs <- t_owsa_ini_pop_costs %>% 
  pivot_longer(1:2, names_to = "group", values_to = "costs") %>% 
  mutate(costs_new = (costs - mod_basecase$total_net_present_cost)) %>% 
  filter(costs_new != 0) %>% 
  mutate(group = ifelse(group == "costs_high", "+25%",
                        ifelse(group == "costs_low", "-25%", group)))


ggplotly(ggplot() +
  geom_bar(data = t_owsa_ini_pop_costs %>% filter(abs(costs_new) >= 500000000),
       aes(x = name, y = costs_new, fill = group, color = group), stat = "identity") + coord_flip() +
  ylab("Difference in Costs from the basecase value") + xlab("") +
  labs(color = "Difference in initial \n population value", fill = "Difference in initial \n population value") +
    scale_x_discrete(labels=c("ini_pop_BN_CHRONIC" = "Chronic pain, no use",
                              "ini_pop_BN_PN" = "Pain free, no use",
                              "ini_pop_BN_ACUTE" = "Acute pain, no use",
                              "ini_pop_BN_CANCER" = "Cancer, no use",
                              "ini_pop_BN_OTHER" = "Other, no use",
                              "ini_pop_BPO_ACUTE" = "Acute pain, Rx use",
                              "ini_pop_BPO_CHRONIC" = "Chronic pain, Rx use",
                              "ini_pop_BPO_CANCER" = "Cancer, Rx use",
                              "ini_pop_BPO_OTHER" = "Other, Rx use",
                              "ini_pop_BPO_MISUSE" = "Rx opioid misuse",
                              "ini_pop_BI_ILLICIT" = "Illicit opioid use")) +
    scale_color_brewer(palette = "Dark2") +
    scale_fill_brewer(palette = "Dark2"))

1.0.0.2 Deaths

Code
t_owsa_ini_pop_deaths <- t_owsa_ini_pop %>% 
  select(deaths_low, deaths_high, name) 
t_owsa_ini_pop_deaths$name <- reorder(t_owsa_ini_pop_deaths$name, t_owsa_ini_pop$deaths_range)

t_owsa_ini_pop_deaths <- t_owsa_ini_pop_deaths %>% 
  pivot_longer(1:2, names_to = "group", values_to = "deaths") %>% 
  mutate(deaths_new = (deaths - mod_basecase$m_M[181, "BO_DEATH"])) %>% 
  filter(deaths_new != 0) %>% 
  mutate(group = ifelse(group == "deaths_high", "+25%",
                        ifelse(group == "deaths_low", "-25%", group)))

ggplotly(ggplot() +
  geom_bar(data = t_owsa_ini_pop_deaths %>% filter(abs(deaths_new) >= 5000),
       aes(x = name, y = deaths_new, fill = group, color = group), stat = "identity") + coord_flip() +
  ylab("Difference in Deaths from the basecase value") + xlab("") +
  labs(color = "Difference in initial \n population value", fill = "Difference in initial \n population value") +
    scale_x_discrete(labels=c("ini_pop_BN_CHRONIC" = "Chronic pain, no use",
                              "ini_pop_BN_PN" = "Pain free, no use",
                              "ini_pop_BN_ACUTE" = "Acute pain, no use",
                              "ini_pop_BN_CANCER" = "Cancer, no use",
                              "ini_pop_BN_OTHER" = "Other, no use",
                              "ini_pop_BPO_ACUTE" = "Acute pain, Rx use",
                              "ini_pop_BPO_CHRONIC" = "Chronic pain, Rx use",
                              "ini_pop_BPO_CANCER" = "Cancer, Rx use",
                              "ini_pop_BPO_OTHER" = "Other, Rx use",
                              "ini_pop_BPO_MISUSE" = "Rx opioid misuse",
                              "ini_pop_BI_ILLICIT" = "Illicit opioid use",
                              "ini_pop_BPO_PALLIATIVE" = "Palliative, Rx use")) +
    scale_color_brewer(palette = "Dark2") +
    scale_fill_brewer(palette = "Dark2"))

1.0.0.3 OD Deaths

Code
t_owsa_ini_pop_oddeaths <- t_owsa_ini_pop %>% 
  select(od_deaths_low, od_deaths_high, name) 
t_owsa_ini_pop_oddeaths$name <- reorder(t_owsa_ini_pop_oddeaths$name, t_owsa_ini_pop$od_deaths_range)

t_owsa_ini_pop_oddeaths <- t_owsa_ini_pop_oddeaths %>% 
  pivot_longer(1:2, names_to = "group", values_to = "oddeaths") %>% 
  mutate(oddeaths_new = (oddeaths - (mod_basecase$m_M[181, "BO_OD_DEATH"] + mod_basecase$extra_od_deaths))) %>% 
  filter(oddeaths_new != 0) %>% 
  mutate(group = ifelse(group == "od_deaths_high", "+25%",
                        ifelse(group == "od_deaths_low", "-25%", group)))

ggplotly(ggplot() +
  geom_bar(data = t_owsa_ini_pop_oddeaths %>% filter(abs(oddeaths_new) >= 100),
       aes(x = name, y = oddeaths_new, fill = group, color = group), stat = "identity") + coord_flip() +
  ylab("Difference in opioid-related overdose deaths from the basecase value") + xlab("") +
  labs(color = "Difference in initial \n population value", fill = "Difference in initial \n population value") +
    scale_x_discrete(labels=c("ini_pop_BN_CHRONIC" = "Chronic pain, no use",
                              "ini_pop_BN_PN" = "Pain free, no use",
                              "ini_pop_BN_ACUTE" = "Acute pain, no use",
                              "ini_pop_BN_CANCER" = "Cancer, no use",
                              "ini_pop_BN_OTHER" = "Other, no use",
                              "ini_pop_BPO_ACUTE" = "Acute pain, Rx use",
                              "ini_pop_BPO_CHRONIC" = "Chronic pain, Rx use",
                              "ini_pop_BPO_CANCER" = "Cancer, Rx use",
                              "ini_pop_BPO_OTHER" = "Other, Rx use",
                              "ini_pop_BPO_MISUSE" = "Rx opioid misuse",
                              "ini_pop_BI_ILLICIT" = "Illicit opioid use",
                              "ini_pop_BPO_PALLIATIVE" = "Palliative, Rx use",
                              "ini_pop_BS_OAT_MAINT" = "OAT maintenance")) +
    scale_color_brewer(palette = "Dark2") +
    scale_fill_brewer(palette = "Dark2"))

2 Transition probabilities

Code
# # gen_arr <- function(m_P_temp) {
# #    arr_P_temp <- array(0, dim=c(n_states, n_states, n_cycles),
# #                dimnames = list(v_state_names, v_state_names, 1:n_cycles))
# #     
# #     p_od_illicit_temp <- m_P_temp["BO_OD_ILLICIT", "BI_ILLICIT"] # 1-rest for transitions from OD
# #     p_illicit_illicit_temp <- m_P_temp["BI_ILLICIT", "BI_ILLICIT"]  # 1-rest for transitions from illicit
# #     
# #     for (t in 1:n_cycles) {
# #       
# #       m_P_this_cycle <- m_P_temp
# #       
# #       if(t %in% year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year %in% c(2015, 2016, 2017)]){
# #         
# #         arr_P_temp[,,t] <- m_P_this_cycle
# #         
# #         } else if (t %in% year_mon_cycle_tbl$cycle[year_mon_cycle_tbl$year %in% c(2018, 2019, 2020)]) {
# #     
# #     m_P_this_cycle["BO_OD_ILLICIT","BO_OD_DEATH"] <- m_P_temp["BO_OD_ILLICIT", "BO_OD_DEATH"] * (1.0055)^(t-36)
# #     m_P_this_cycle["BO_OD_ILLICIT", "BI_ILLICIT"] <- 1-(rowSums(m_P_this_cycle)["BO_OD_ILLICIT"] - p_od_illicit_temp)
# #     
# #     m_P_this_cycle["BI_ILLICIT", "BO_OD_ILLICIT"] <- m_P_temp["BI_ILLICIT", "BO_OD_ILLICIT"] * (1.0055)^(t-36)
# #     m_P_this_cycle["BI_ILLICIT", "BI_ILLICIT"] <- 1-(rowSums(m_P_this_cycle)["BI_ILLICIT"] - p_illicit_illicit_temp)
# #   
# #   } else {
# #     
# #     m_P_this_cycle["BO_OD_ILLICIT","BO_OD_DEATH"] <- m_P_temp["BO_OD_ILLICIT", "BO_OD_DEATH"] * (1.0055)^(36)
# #     m_P_this_cycle["BO_OD_ILLICIT", "BI_ILLICIT"] <- 1-(rowSums(m_P_this_cycle)["BO_OD_ILLICIT"] - p_od_illicit_temp)
# #     
# #     m_P_this_cycle["BI_ILLICIT", "BO_OD_ILLICIT"] <- m_P_temp["BI_ILLICIT", "BO_OD_ILLICIT"] * (1.0055)^(36)
# #     m_P_this_cycle["BI_ILLICIT", "BI_ILLICIT"] <- 1-(rowSums(m_P_this_cycle)["BI_ILLICIT"] - p_illicit_illicit_temp)
# #      
# #   }
#    
# #   arr_P_temp[,,t] <- m_P_this_cycle
# #     }
# #     return(arr_P_temp)
# # }

# # 
# # trans_prob_tbl_new <- trans_prob_tbl_new %>%
#   full_join(., trans_prob_tbl %>%
# #               select(var_from, var_to, lb, ub), by = c("var_from", "var_to")) %>%
# #   replace(is.na(.), 0)
Code
# trans_prob_tbl_ows <- cbind.data.frame(rep(v_state_names, each = n_states), 
#                                          rep(v_state_names, n_states))
# 
# names(trans_prob_tbl_ows) <- c("var_from", "var_to")
# names(trans_prob_tbl) <- c("var_from", "var_to",
#                  "basevalue", "lb", "ub", "range")
# 
# trans_prob_tbl_ows <- trans_prob_tbl_ows %>% 
#     full_join(., trans_prob_tbl %>% 
#                 select(var_from, var_to, basevalue, lb, ub), by = c("var_from", "var_to")) %>% 
#     replace(is.na(.), 0)
# 
# v_var1 <- unique(trans_prob_tbl_ows$var_from[!(trans_prob_tbl_ows$var_from %in% c("BO_OD_DEATH", "BO_DEATH"))])
# v_var2 <- trans_prob_tbl_ows$var_to[trans_prob_tbl_ows$basevalue == 999]
# 
# 
# dsa_fun <- function(var1, var2) {
#   
#   m_temp <- matrix(0,
#                    nrow = nrow(trans_prob_tbl_ows[trans_prob_tbl_ows$var_from == var1, ]),
#                    ncol = 8 + 4,
#                    dimnames = list(
#                      rows = paste0("p_", trans_prob_tbl_ows$var_to[trans_prob_tbl_ows$var_from == var1]),
#                      cols = c("input_bc", "input_low", "input_high",
#                               "costs_low", "costs_high", "cost_bc",
#                               "deaths_low", "deaths_high", "deaths_bc",
#                               "od_deaths_low", "od_deaths_high", "od_deaths_bc")
#                          ))
#   
#   for(state_to in unique(trans_prob_tbl_ows$var_to)){
#     
#     t_p_temp_full <- trans_prob_tbl_ows %>% 
#       mutate(basevalue = ifelse(var_to == state_to & 
#                                   var_from == var1 &
#                                   var_to != var2, lb, basevalue))
# 
#     m_P_temp <- matrix(t_p_temp_full$basevalue, byrow = T,
#                 nrow = n_states, ncol = n_states,
#                 dimnames = list(v_state_names, v_state_names))
#     
#     for (i in 1:n_states){
#       for (j in 1:n_states){
#         m_P_temp[i, j] <- ifelse(m_P_temp[i, j] == 999, 1 - (rowSums(m_P_temp)[i] - 999), m_P_temp[i, j])
#       }
#     }
#     
#     arr_P_temp <- gen_trans_prob_arr(m_P_temp)
#     
#     mod <- markov_mod(array_P = arr_P_temp, mat_P = NA)
#     
#     deaths_lb <- mod$m_M[181, "BO_DEATH"]
#     od_deaths_lb <- mod$m_M[181, "BO_OD_DEATH"] + mod$extra_od_deaths
#     costs_lb <- mod$total_net_present_cost
#     
#     t_p_temp_full <- trans_prob_tbl_ows %>% 
#       mutate(basevalue = ifelse(var_to == state_to & 
#                                   var_from == var1 &
#                                   var_to != var2, ub, basevalue))
#     
#     m_P_temp <- matrix(t_p_temp_full$basevalue, byrow = T,
#               nrow = n_states, ncol = n_states,
#               dimnames = list(v_state_names, v_state_names))
#     
#     
#     for (i in 1:n_states){
#       for (j in 1:n_states){
#         m_P_temp[i, j] <- ifelse(m_P_temp[i, j] == 999, 1 - (rowSums(m_P_temp)[i] - 999), m_P_temp[i, j])
#       }
#     }
#     
#     arr_P_temp <- gen_trans_prob_arr(m_P_temp)
#     
#     mod <- markov_mod(array_P = arr_P_temp, mat_P = NA)
#     
#     
#     deaths_ub <- mod$m_M[181, "BO_DEATH"]
#     od_deaths_ub <- mod$m_M[181, "BO_OD_DEATH"] + mod$extra_od_deaths
#     costs_ub <- mod$total_net_present_cost
#     
#     
#     t_p_temp_full <- trans_prob_tbl_ows
#     
#     m_P_temp <- matrix(t_p_temp_full$basevalue, byrow = T,
#               nrow = n_states, ncol = n_states,
#               dimnames = list(v_state_names, v_state_names))
#     
#     
#     for (i in 1:n_states){
#       for (j in 1:n_states){
#         m_P_temp[i, j] <- ifelse(m_P_temp[i, j] == 999, 1 - (rowSums(m_P_temp)[i] - 999), m_P_temp[i, j])
#       }
#     }
#     
#     arr_P_temp <- gen_trans_prob_arr(m_P_temp)
#     
#     mod <- markov_mod(array_P = arr_P_temp, mat_P = NA)
#     
#     
#     deaths_bc_ows <- mod$m_M[181, "BO_DEATH"]
#     od_deaths_bc_ows <- mod$m_M[181, "BO_OD_DEATH"] + mod$extra_od_deaths
#     costs_bc_ows <- mod$total_net_present_cost
#     
#     
#     
#     
#     m_temp[paste0("p_", state_to), ] <- c(trans_prob_tbl_ows$basevalue[trans_prob_tbl_ows$var_from == var1 &
#                                                                   trans_prob_tbl_ows$var_to == state_to],
#                                           trans_prob_tbl_ows$lb[trans_prob_tbl_ows$var_from == var1 &
#                                                                   trans_prob_tbl_ows$var_to == state_to],
#                                           trans_prob_tbl_ows$ub[trans_prob_tbl_ows$var_from == var1 &
#                                                                   trans_prob_tbl_ows$var_to == state_to],
#                                           costs_lb, costs_ub,costs_bc_ows,
#                                           deaths_lb, deaths_ub,deaths_bc_ows,
#                                           od_deaths_lb, od_deaths_ub, od_deaths_bc_ows)
#   }
#   
#   return(as_tibble(m_temp) %>% 
#   mutate(group = var1,
#          cost_range = abs(costs_high - costs_low),
#          deaths_range = abs(deaths_high - deaths_low),
#          od_deaths_range = abs(od_deaths_high - od_deaths_low),
#          name = rownames(m_temp)))
#   
# }
Code
# # basecase for ows
# 
# # m_P_bc <- matrix(trans_prob_tbl_new$basevalue, byrow = T,
# #               nrow = n_states, ncol = n_states,
# #               dimnames = list(v_state_names, v_state_names))
# #     
# #     for (i in 1:n_states){
# #       for (j in 1:n_states){
# #         m_P_bc[i, j] <- ifelse(m_P_bc[i, j] == 999, 1 - (rowSums(m_P_bc)[i] - 999), m_P_bc[i, j])
# #       }
# #     }
#  
# # 
# # 
# #     
# # modbc <- markov_mod(array_P = NA, mat_P = m_P_bc)
# # deaths_bc <- modbc$m_M[181, "BO_DEATH"]
# # od_deaths_bc <- modbc$m_M[181, "BO_OD_DEATH"] + modbc$extra_od_deaths
# # costs_bc <- modbc$total_net_present_cost
# # 
# # 
# # dsa_fun <- function(var1, var2) {
# #   
# #   m_temp <- matrix(0,
# #                    nrow = nrow(trans_prob_tbl_new[trans_prob_tbl_new$var_from == var1, ]),
# #                    ncol = 8 + 4,
# #                    dimnames = list(
# #                      rows = paste0("p_", trans_prob_tbl_new$var_to[trans_prob_tbl_new$var_from == var1]),
# #                      cols = c("input_bc", "input_low", "input_high",
# #                               "costs_low", "costs_high", "cost_bc",
# #                               "deaths_low", "deaths_high", "deaths_bc",
# #                               "od_deaths_low", "od_deaths_high", "od_deaths_bc")
# #                          ))
# #   
# #   for(state_to in unique(trans_prob_tbl_new$var_to)){
# #     
# #     t_p_temp_full <- trans_prob_tbl_new %>% 
# #       mutate(basevalue = ifelse(var_to == state_to & 
# #                                   var_from == var1 &
# #                                   var_to != var2, lb, basevalue))
# #     
# #     m_P_temp <- matrix(t_p_temp_full$basevalue, byrow = T,
# #               nrow = n_states, ncol = n_states,
# #               dimnames = list(v_state_names, v_state_names))
# #     
# #     for (i in 1:n_states){
# #       for (j in 1:n_states){
# #         m_P_temp[i, j] <- ifelse(m_P_temp[i, j] == 999, 1 - (rowSums(m_P_temp)[i] - 999), m_P_temp[i, j])
# #       }
# #       }
# #     
# #     mod <- markov_mod(array_P = NA, mat_P = m_P_temp)
# #     
# #     deaths_lb <- mod$m_M[181, "BO_DEATH"]
# #     od_deaths_lb <- mod$m_M[181, "BO_OD_DEATH"] + mod$extra_od_deaths
# #     costs_lb <- mod$total_net_present_cost
# #     
# #     
# #      t_p_temp_full <- trans_prob_tbl_new %>% 
# #       mutate(basevalue = ifelse(var_to == state_to & 
# #                                   var_from == var1 &
# #                                   var_to != var2, ub, basevalue))
# #     
# #     m_P_temp <- matrix(t_p_temp_full$basevalue, byrow = T,
# #               nrow = n_states, ncol = n_states,
# #               dimnames = list(v_state_names, v_state_names))
# #     
# #     
# #     for (i in 1:n_states){
# #       for (j in 1:n_states){
# #         m_P_temp[i, j] <- ifelse(m_P_temp[i, j] == 999, 1 - (rowSums(m_P_temp)[i] - 999), m_P_temp[i, j])
# #       }
# #       }
# #     
# #     mod <- markov_mod(array_P = NA, mat_P = m_P_temp)
# #     
# #     
# #     deaths_ub <- mod$m_M[181, "BO_DEATH"]
# #     od_deaths_ub <- mod$m_M[181, "BO_OD_DEATH"] + mod$extra_od_deaths
# #     costs_ub <- mod$total_net_present_cost
# #     
# #     
# #     t_p_temp_full <- trans_prob_tbl_new
# #     
# #     m_P_temp <- matrix(t_p_temp_full$basevalue, byrow = T,
# #               nrow = n_states, ncol = n_states,
# #               dimnames = list(v_state_names, v_state_names))
# #     
# #     for (i in 1:n_states){
# #       for (j in 1:n_states){
# #         m_P_temp[i, j] <- ifelse(m_P_temp[i, j] == 999, 1 - (rowSums(m_P_temp)[i] - 999), m_P_temp[i, j])
# #       }
# #       }
# #     
# #     mod <- markov_mod(array_P = NA, mat_P = m_P_temp)
# #     
# #     deaths_bc_tmp <- mod$m_M[181, "BO_DEATH"]
# #     od_deaths_bc_tmp <- mod$m_M[181, "BO_OD_DEATH"] + mod$extra_od_deaths
# #     costs_bc_tmp <- mod$total_net_present_cost
# #     
# #     
# #     m_temp[paste0("p_", state_to), ] <- c(trans_prob_tbl_new$basevalue[trans_prob_tbl_new$var_from == var1 &
# #                                                                   trans_prob_tbl_new$var_to == state_to],
# #                                           trans_prob_tbl_new$lb[trans_prob_tbl_new$var_from == var1 &
# #                                                                   trans_prob_tbl_new$var_to == state_to],
# #                                           trans_prob_tbl_new$ub[trans_prob_tbl_new$var_from == var1 &
# #                                                                   trans_prob_tbl_new$var_to == state_to],
# #                                           costs_lb, costs_ub,costs_bc_tmp,
# #                                           deaths_lb, deaths_ub,deaths_bc_tmp,
# #                                           od_deaths_lb, od_deaths_ub, od_deaths_bc_tmp)
# #   }
# #   
# #   return(as_tibble(m_temp) %>% 
# #   mutate(group = var1,
# #          cost_range = abs(costs_high - costs_low),
# #          deaths_range = abs(deaths_high - deaths_low),
# #          od_deaths_range = abs(od_deaths_high - od_deaths_low),
# #          name = rownames(m_temp)))
# #   
# # }
# #
Code
# ows_tbl_prob <- data.frame()
# 
# for (i in 1:length(v_var1)){
#   ows_tbl_prob <- rbind.data.frame(ows_tbl_prob, dsa_fun(v_var1[i], v_var2[i]))
# }
# 
# write.csv(ows_tbl_prob, file = "../01_data/ows_tbl_prob.csv")

2.0.0.1 Costs

Code
ows_tbl_prob_costs <- ows_tbl_prob %>% 
  select(costs_low, costs_high, name, group) %>%
  mutate(name = paste(name, group, sep = "_")) %>% 
  rename(grp = group)

ows_tbl_prob_costs$name <- reorder(ows_tbl_prob_costs$name, ows_tbl_prob$cost_range)

ows_tbl_prob_costs <- ows_tbl_prob_costs %>% 
  pivot_longer(1:2, names_to = "group", values_to = "costs") %>% 
  mutate(costs_new = (costs - as.numeric(mod_basecase$total_net_present_cost)),
         costs_new_per = (costs_new/as.numeric(mod_basecase$total_net_present_cost)) * 100) %>% 
  filter(costs_new != 0) %>% 
  mutate(group = ifelse(group == "costs_high", "+25%",
                        ifelse(group == "costs_low", "-25%", group)))


name_cost_1 <- ows_tbl_prob_costs %>% filter(abs(costs_new_per) >= 0.5) %>% select(name)
name_cost_2 <- ows_tbl_prob_costs %>% filter(abs(costs_new_per) >= 1.5) %>% select(name)

ggplotly(ggplot() +
  geom_bar(data = ows_tbl_prob_costs %>% filter(name %in% name_cost_1$name),
       aes(x = name, y = costs_new, fill = group, color = group), stat = "identity") + coord_flip() +
  ylab("Difference in Costs from the basecase value") + xlab("") +
  labs(color = "Difference in transition \n probability value", fill = "Difference in transition \n probability value") +
    scale_x_discrete(labels=c("p_BN_ACUTE_BN_PN" = "Pain free, no use to acute pain, no use",
                              "p_BN_CANCER_BN_PN" = "Pain free, no use to cancer, no use",
                              "p_BPO_ACUTE_BN_PN" = "Pain free, no use to acute pain, Rx use",
                              "p_BPO_MISUSE_BN_PN" = "Pain free, no use to Rx opioid misuse",
                              "p_BN_CHRONIC_BN_ACUTE" = "Acute pain, no use to chronic pain, no use",
                              "p_BN_CANCER_BN_CHRONIC" = "Chronic pain, no use to cancer, no use",
                              "p_BPO_MISUSE_BN_CHRONIC" = "Chronic pain, no use to Rx opioid misuse",
                              "p_BPO_OTHER_BN_OTHER" = "Other, no use to other, Rx use",
                              "p_BN_CHRONIC_BPO_ACUTE" = "Acute pain, Rx use to chronic pain, no use",
                              "p_BPO_MISUSE_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid misuse",
                              "p_BO_OD_RX_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid overdose",
                              "p_BO_DEATH_BPO_CANCER" = "Cancer, Rx use to death",
                              "p_BN_OTHER_BPO_OTHER" = "Other, Rx use to other, no use",
                              "p_BPO_MISUSE_BPO_OTHER" = "Other, Rx use to Rx opioid misuse",
                              "p_BN_PN_BPO_MISUSE" = "Rx opioid misuse to pain free, no use",
                              "p_BI_ILLICIT_BPO_MISUSE" = "Rx opioid misuse to illicit opioid use",
                              "p_BS_DETOX_BI_ILLICIT" = "Illicit opioid use to detox/withdrawal man.",
                              "p_BS_OAT_INI_BS_DETOX" = "Detox/withdrawal man. to OAT initiation",
                              "p_BS_OAT_MAINT_BS_OAT_INI" = "OAT initiation to OAT maintenance",
                              "p_BS_OAT_MAINT_BS_OAT_MAINT" = "OAT maintenance stay",
                              "p_BO_OD_DEATH_BO_OD_RX" = "Rx opioid overdose to death",
                              "p_BN_PALLIATIVE_BN_PN" = "Pain Free, no use, to Palliative, no use",
                              "p_BO_DEATH_BN_PN" = "Pain Free, no use, to Death",
                              "p_BN_PALLIATIVE_BN_CHRONIC" = "Chronic pain, no use to palliative, no use",
                              "p_BPO_CHRONIC_BN_CHRONIC" = "Chronic pain, no use to chronic pain, Rx use",  
                              "p_BO_DEATH_BN_CHRONIC" = "Chronic pain, no use to death",
                              "p_BO_DEATH_BN_OTHER" = "Other no use to death",
                              "p_BN_CHRONIC_BPO_CHRONIC" = "Chronic pain, Rx use to chronic pain, no use",
                              "p_BO_DEATH_BPO_CHRONIC" = "Chronic pain, Rx use to death", 
                              "p_BPO_MISUSE_BPO_CANCER" = "Cancer, Rx use to Rx opioid misuse",
                              "p_BO_OD_RX_BPO_CANCER" = "Cancer, Rx use to Rx opioid overdose",
                              "p_BO_DEATH_BPO_MISUSE" = "Rx opioid misuse to death",
                              "p_BR_OAT_MAINT_BR_OAT_MAINT" = "R: OAT maintenance stay")) +
    scale_color_brewer(palette = "Dark2") +
    scale_fill_brewer(palette = "Dark2"))
Code
ggplotly(ggplot() +
  geom_bar(data = ows_tbl_prob_costs %>% filter(name %in% name_cost_2$name),
       aes(x = name, y = costs_new, fill = group, color = group), stat = "identity") + coord_flip() +
  ylab("Difference in Costs from the basecase value") + xlab("") +
  labs(color = "Difference in transition \n probability value", fill = "Difference in transition \n probability value") +
    scale_x_discrete(labels=c("p_BN_ACUTE_BN_PN" = "Pain free, no use to acute pain, no use",
                              "p_BN_CANCER_BN_PN" = "Pain free, no use to cancer, no use",
                              "p_BPO_ACUTE_BN_PN" = "Pain free, no use to acute pain, Rx use",
                              "p_BPO_MISUSE_BN_PN" = "Pain free, no use to Rx opioid misuse",
                              "p_BN_CHRONIC_BN_ACUTE" = "Acute pain, no use to chronic pain, no use",
                              "p_BN_CANCER_BN_CHRONIC" = "Chronic pain, no use to cancer, no use",
                              "p_BPO_MISUSE_BN_CHRONIC" = "Chronic pain, no use to Rx opioid misuse",
                              "p_BPO_OTHER_BN_OTHER" = "Other, no use to other, Rx use",
                              "p_BN_CHRONIC_BPO_ACUTE" = "Acute pain, Rx use to chronic pain, no use",
                              "p_BPO_MISUSE_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid misuse",
                              "p_BO_OD_RX_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid overdose",
                              "p_BO_DEATH_BPO_CANCER" = "Cancer, Rx use to death",
                              "p_BN_OTHER_BPO_OTHER" = "Other, Rx use to other, no use",
                              "p_BPO_MISUSE_BPO_OTHER" = "Other, Rx use to Rx opioid misuse",
                              "p_BN_PN_BPO_MISUSE" = "Rx opioid misuse to pain free, no use",
                              "p_BI_ILLICIT_BPO_MISUSE" = "Rx opioid misuse to illicit opioid use",
                              "p_BS_DETOX_BI_ILLICIT" = "Illicit opioid use to detox/withdrawal man.",
                              "p_BS_OAT_INI_BS_DETOX" = "Detox/withdrawal man. to OAT initiation",
                              "p_BS_OAT_MAINT_BS_OAT_INI" = "OAT initiation to OAT maintenance",
                              "p_BS_OAT_MAINT_BS_OAT_MAINT" = "OAT maintenance stay",
                              "p_BO_OD_DEATH_BO_OD_RX" = "Rx opioid overdose to death",
                              "p_BN_PALLIATIVE_BN_PN" = "Pain Free, no use, to Palliative, no use",
                              "p_BO_DEATH_BN_PN" = "Pain Free, no use, to Death",
                              "p_BN_PALLIATIVE_BN_CHRONIC" = "Chronic pain, no use to palliative, no use",
                              "p_BPO_CHRONIC_BN_CHRONIC" = "Chronic pain, no use to chronic pain, Rx use",  
                              "p_BO_DEATH_BN_CHRONIC" = "Chronic pain, no use to death",
                              "p_BO_DEATH_BN_OTHER" = "Other no use to death",
                              "p_BN_CHRONIC_BPO_CHRONIC" = "Chronic pain, Rx use to chronic pain, no use",
                              "p_BO_DEATH_BPO_CHRONIC" = "Chronic pain, Rx use to death", 
                              "p_BPO_MISUSE_BPO_CANCER" = "Cancer, Rx use to Rx opioid misuse",
                              "p_BO_OD_RX_BPO_CANCER" = "Cancer, Rx use to Rx opioid overdose",
                              "p_BO_DEATH_BPO_MISUSE" = "Rx opioid misuse to death",
                              "p_BR_OAT_MAINT_BR_OAT_MAINT" = "R: OAT maintenance stay")) +
    scale_color_brewer(palette = "Dark2") +
    scale_fill_brewer(palette = "Dark2"))

2.0.0.2 Deaths

Code
ows_tbl_prob_deaths <- ows_tbl_prob %>% 
  select(deaths_low, deaths_high, name, group) %>%
  mutate(name = paste(name, group, sep = "_")) %>% 
  rename(grp = group)

ows_tbl_prob_deaths$name <- reorder(ows_tbl_prob_deaths$name, ows_tbl_prob$deaths_range)

ows_tbl_prob_deaths <- ows_tbl_prob_deaths %>% 
  pivot_longer(1:2, names_to = "group", values_to = "deaths") %>% 
  mutate(deaths_new = (deaths - as.numeric(mod_basecase$m_M[181, "BO_DEATH"])),
         deaths_new_per = (deaths_new/as.numeric(mod_basecase$m_M[181, "BO_DEATH"])) * 100) %>% 
  filter(deaths_new != 0) %>% 
  mutate(group = ifelse(group == "deaths_high", "+25%",
                        ifelse(group == "deaths_low", "-25%", group)))


name_death_1 <- ows_tbl_prob_deaths %>% filter(abs(deaths_new_per) >= 0.5) %>% select(name)

ggplotly(ggplot() +
  geom_bar(data = ows_tbl_prob_deaths %>% filter(name %in% name_death_1$name),
       aes(x = name, y = deaths_new, fill = group, color = group), stat = "identity") + coord_flip() +
  ylab("Difference in deaths from the basecase value") + xlab("") +
  labs(color = "Difference in transition \n probability value", fill = "Difference in transition \n probability value") +
    scale_x_discrete(labels=c("p_BN_ACUTE_BN_PN" = "Pain free, no use to acute pain, no use",
                              "p_BN_CANCER_BN_PN" = "Pain free, no use to cancer, no use",
                              "p_BPO_ACUTE_BN_PN" = "Pain free, no use to acute pain, Rx use",
                              "p_BPO_MISUSE_BN_PN" = "Pain free, no use to Rx opioid misuse",
                              "p_BN_CHRONIC_BN_ACUTE" = "Acute pain, no use to chronic pain, no use",
                              "p_BN_CANCER_BN_CHRONIC" = "Chronic pain, no use to cancer, no use",
                              "p_BPO_MISUSE_BN_CHRONIC" = "Chronic pain, no use to Rx opioid misuse",
                              "p_BPO_OTHER_BN_OTHER" = "Other, no use to other, Rx use",
                              "p_BN_CHRONIC_BPO_ACUTE" = "Acute pain, Rx use to chronic pain, no use",
                              "p_BPO_MISUSE_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid misuse",
                              "p_BO_OD_RX_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid overdose",
                              "p_BO_DEATH_BPO_CANCER" = "Cancer, Rx use to death",
                              "p_BN_OTHER_BPO_OTHER" = "Other, Rx use to other, no use",
                              "p_BPO_MISUSE_BPO_OTHER" = "Other, Rx use to Rx opioid misuse",
                              "p_BN_PN_BPO_MISUSE" = "Rx opioid misuse to pain free, no use",
                              "p_BI_ILLICIT_BPO_MISUSE" = "Rx opioid misuse to illicit opioid use",
                              "p_BS_DETOX_BI_ILLICIT" = "Illicit opioid use to detox/withdrawal man.",
                              "p_BS_OAT_INI_BS_DETOX" = "Detox/withdrawal man. to OAT initiation",
                              "p_BS_OAT_MAINT_BS_OAT_INI" = "OAT initiation to OAT maintenance",
                              "p_BS_OAT_MAINT_BS_OAT_MAINT" = "OAT maintenance stay",
                              "p_BO_OD_DEATH_BO_OD_RX" = "Rx opioid overdose to death",
                              "p_BN_PALLIATIVE_BN_PN" = "Pain Free, no use, to Palliative, no use",
                              "p_BO_DEATH_BN_PN" = "Pain Free, no use, to Death",
                              "p_BN_PALLIATIVE_BN_CHRONIC" = "Chronic pain, no use to palliative, no use",
                              "p_BPO_CHRONIC_BN_CHRONIC" = "Chronic pain, no use to chronic pain, Rx use",  
                              "p_BO_DEATH_BN_CHRONIC" = "Chronic pain, no use to death",
                              "p_BO_DEATH_BN_OTHER" = "Other no use to death",
                              "p_BN_CHRONIC_BPO_CHRONIC" = "Chronic pain, Rx use to chronic pain, no use",
                              "p_BO_DEATH_BPO_CHRONIC" = "Chronic pain, Rx use to death", 
                              "p_BPO_MISUSE_BPO_CANCER" = "Cancer, Rx use to Rx opioid misuse",
                              "p_BO_OD_RX_BPO_CANCER" = "Cancer, Rx use to Rx opioid overdose",
                              "p_BO_DEATH_BPO_MISUSE" = "Rx opioid misuse to death",
                              "p_BR_OAT_MAINT_BR_OAT_MAINT" = "R: OAT maintenance stay")) +
    scale_color_brewer(palette = "Dark2") +
    scale_fill_brewer(palette = "Dark2"))

2.0.0.3 OD-deaths

Code
oddeaths_bc <- as.numeric(mod_basecase$m_M[181, "BO_OD_DEATH"] + mod_basecase$extra_od_deaths) 
ows_tbl_prob_oddeaths <- ows_tbl_prob %>% 
  select(od_deaths_low, od_deaths_high, name, group) %>%
  mutate(name = paste(name, group, sep = "_")) %>% 
  rename(grp = group)

ows_tbl_prob_oddeaths$name <- reorder(ows_tbl_prob_oddeaths$name, ows_tbl_prob$od_deaths_range)

ows_tbl_prob_oddeaths <- ows_tbl_prob_oddeaths %>% 
  pivot_longer(1:2, names_to = "group", values_to = "oddeaths") %>% 
  mutate(oddeaths_new = (oddeaths - oddeaths_bc),
         oddeaths_new_per = (oddeaths_new/oddeaths_bc) * 100) %>% 
  filter(oddeaths_new != 0) %>% 
  mutate(group = ifelse(group == "od_deaths_high", "+25%",
                        ifelse(group == "od_deaths_low", "-25%", group)))


name_oddeath_1 <- ows_tbl_prob_oddeaths %>% filter(abs(oddeaths_new_per) >= 5) %>% select(name)

ggplotly(ggplot() +
  geom_bar(data = ows_tbl_prob_oddeaths %>% filter(name %in% name_oddeath_1$name),
       aes(x = name, y = oddeaths_new, fill = group, color = group), stat = "identity") + coord_flip() +
  ylab("Difference in opioid-related overdose death from the basecase value") + xlab("") +
  labs(color = "Difference in transition \n probability value", fill = "Difference in transition \n probability value") +
    scale_x_discrete(labels=c("p_BN_ACUTE_BN_PN" = "Pain free, no use to acute pain, no use",
                              "p_BN_CANCER_BN_PN" = "Pain free, no use to cancer, no use",
                              "p_BPO_ACUTE_BN_PN" = "Pain free, no use to acute pain, Rx use",
                              "p_BPO_MISUSE_BN_PN" = "Pain free, no use to Rx opioid misuse",
                              "p_BN_CHRONIC_BN_ACUTE" = "Acute pain, no use to chronic pain, no use",
                              "p_BN_CANCER_BN_CHRONIC" = "Chronic pain, no use to cancer, no use",
                              "p_BPO_MISUSE_BN_CHRONIC" = "Chronic pain, no use to Rx opioid misuse",
                              "p_BPO_OTHER_BN_OTHER" = "Other, no use to other, Rx use",
                              "p_BN_CHRONIC_BPO_ACUTE" = "Acute pain, Rx use to chronic pain, no use",
                              "p_BPO_MISUSE_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid misuse",
                              "p_BO_OD_RX_BPO_CHRONIC" = "Chronic pain, Rx use to Rx opioid overdose",
                              "p_BO_DEATH_BPO_CANCER" = "Cancer, Rx use to death",
                              "p_BN_OTHER_BPO_OTHER" = "Other, Rx use to other, no use",
                              "p_BPO_MISUSE_BPO_OTHER" = "Other, Rx use to Rx opioid misuse",
                              "p_BN_PN_BPO_MISUSE" = "Rx opioid misuse to pain free, no use",
                              "p_BI_ILLICIT_BPO_MISUSE" = "Rx opioid misuse to illicit opioid use",
                              "p_BS_DETOX_BI_ILLICIT" = "Illicit opioid use to detox/withdrawal man.",
                              "p_BS_OAT_INI_BS_DETOX" = "Detox/withdrawal man. to OAT initiation",
                              "p_BS_OAT_MAINT_BS_OAT_INI" = "OAT initiation to OAT maintenance",
                              "p_BS_OAT_MAINT_BS_OAT_MAINT" = "OAT maintenance stay",
                              "p_BO_OD_DEATH_BO_OD_RX" = "Rx opioid overdose to death",
                              "p_BN_PALLIATIVE_BN_PN" = "Pain Free, no use, to Palliative, no use",
                              "p_BO_DEATH_BN_PN" = "Pain Free, no use, to Death",
                              "p_BN_PALLIATIVE_BN_CHRONIC" = "Chronic pain, no use to palliative, no use",
                              "p_BPO_CHRONIC_BN_CHRONIC" = "Chronic pain, no use to chronic pain, Rx use",  
                              "p_BO_DEATH_BN_CHRONIC" = "Chronic pain, no use to death",
                              "p_BO_DEATH_BN_OTHER" = "Other no use to death",
                              "p_BN_CHRONIC_BPO_CHRONIC" = "Chronic pain, Rx use to chronic pain, no use",
                              "p_BO_DEATH_BPO_CHRONIC" = "Chronic pain, Rx use to death", 
                              "p_BPO_MISUSE_BPO_CANCER" = "Cancer, Rx use to Rx opioid misuse",
                              "p_BO_OD_RX_BPO_CANCER" = "Cancer, Rx use to Rx opioid overdose",
                              "p_BO_DEATH_BPO_MISUSE" = "Rx opioid misuse to death",
                              "p_BR_OAT_MAINT_BR_OAT_MAINT" = "R: OAT maintenance stay",
                              "p_BO_OD_DEATH_BO_OD_ILLICIT" = "Illicit opioid overdose to death",
                              "p_BO_OD_ILLICIT_BI_ILLICIT" = "Illicit opioid use to Illicit opioid overdose",
                              "p_BO_OD_RX_BPO_MISUSE" = "Rx opioid misuse to Rx opioid overdose",
                              "p_BR_OD_ILLICIT_BR_ILLICIT" = "R: illicit opioid use to R: illicit opioid overdose",
                              "p_BO_OD_DEATH_BR_OD_ILLICIT" = "R: Illicit opioid overdose to death",
                              "p_BO_OD_ILLICIT_BS_DETOX" = "Detox/withdrawal man. to Illicit opioid overdose",
                              "p_BI_ILLICIT_BS_DETOX" = "Detox/withdrawal man. to illicit opioid use")) +
    scale_color_brewer(palette = "Dark2") +
    scale_fill_brewer(palette = "Dark2"))